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Abstract 

(3 ■ This paper extends the Chandrasekhar-type recursions due to Morf, Sidhu, and Kailath "Some new 

(N' 

>• ! algorithms for recursive estimation in constant, linear, discrete-time systems, IEEE Trans. Autom. 

O 



Control 19 (1974) 315-323" to the case of periodic time- varying state-space models. We show 
that the S'-lagged increments of the one-step prediction error covariance satisfy certain recursions 
from which we derive some algorithms for linear least squares estimation for periodic state-space 
models. The proposed recursions may have potential computational advantages over the Kalman 
Filter and, in particular, the periodic Riccati difference equation. 



r^ ' Keywords: Periodic state-space models, Chandrasekhar-type recursions, Kalman Filter, periodic 

in 

00 ■ Riccati difference equation. 

t^^ ' Introduction 

O 

Morf et al (1974) proposed recursions that substitute the Kalman Filter for linear least squares 
3 . estimation of discrete-time time-invariant state space models, with a simpler computational complex- 

ity. The new algorithms have been called Chandrasekhar-type recursions because they are analog 
to certain differential equations encountered in continuous -time problems (Kailath, 1973). Since 
there, a considerable attention has been paid in the three recent decades to the Chandrasekhar-type 
recursions (see e.g. Friedlander et al, 1978; Morf and Kailath; 1975; Houacine and Demonent, 
1986; Houacine, 1991; Nakamori et al, 2004; Nakamori, 2007). At present, there exist several 
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useful applications of the Chandrasekhar filter in improving computational aspects related to the 
building of linear time-invariant models. We mention non exhaustively the likelihood evaluation 
(see Pearlman, 1980; Melard, 1984; Kohn and Ansley, 1985 for autoregressive moving average, 
ARM A, models and Shea, 1989 for vector ARM A models), the calculation of the exact Fisher 
information matrix (see Melard and Klein (1994) for the ARM A case and Klein et al (1998) 
for general dynamic time-invariant models), and the development of fast variants of the recursive 
least squares algorithm (Houacine, 1991; Sayed and Kailath, 1994; Nakamori et al 2004). As is 
well known, the Chandrasekhar equations are restricted to the case of time-invariant state-space 
models because of their particular time invariance structure and it seems that there is no results 
tied to the class of all nonstationarity, except in very special cases (Sayed and Kailath, 1994). 
A particular class of nonstationarity whose importance has no need to be proven is the one of 
periodic linear models. Important progress has been made recently in the building and analysis 
of periodic ARM A (PARMA) and periodic state-space characterizations. The objective was to 
develop extensions of similar methods for standard time-invariant models to their periodic coun- 
terparts, without transforming periodic systems to their corresponding multivariate time-invariant 
representations in order to simplify the computational burden. Despite the current abundance of 
computational methods for periodic state-space models (see e.g. Lund and Basawa, 2000; Varga 
and Van Dooren, 2001; Gautier, 2005; Bentarzi and Aknouche, 2005; Aknouche, 2007; Aknouche 
and Hamdi, 2007 Aknouche et al, 2007 and the references therein) it seems that there is no results 
concerning extensions of the Chandrasekhar recursions to the periodic case. This paper proposes 
some algorithms for linear least squares estimation of periodic state-space models. Our methods 
extend the Chandrasekhar algorithms proposed by Morf et al (1974) to the periodic time- varying 
case and retain their desirable features. As a result, the periodic Chandrasekhar recursions are used 
through the innovation approach to efficiently evaluate the likelihood of periodic ARM A models. 
The rest of this paper is organized as follows. Section 1 briefly recalls some preliminary defini- 



tions and facts about periodic state-space models and their corresponding Kalman Filter. In Section 
2 we develop some Chandrasekhar-type algorithms that substitute the Kalman filter for periodic 
state-space models. The initialization problem will be studied in Section 3. 



I. Preliminary definitions and notations 



Consider the following linear periodic state-space model 



xt+i = FtXf + Gtet 



t e Z, 



(1) 



where {xi}, {yt}, {et} and {et} are random processes of dimensions r x 1, m x 1, m x 1, and 
d X I respectively, with 

E{et) = E{et)=0 



E {ett[^h) = Sh,oQt and < 



E {etK-k) = 
E (etyU) = 



\/t,heZ 



Wk>0 



(S stands for the Kronecker function). The nonrandom matrices Ft, Gt, H[, Qt, Rt, and Wt are 
periodic in time with period 5". To simplify the exposition we suppose without loss of generality 
that 

E (ete'i) = 0, Vt, / G Z. 

Let Xi and y^ be the linear least squares forecasts of x^ and y^, respectively, based on xi, X2..., Xi_i. 
Then as is well known, x^ and y^ may be uniquely obtained from the Kalman filter (Kalman, 1960) 
which is given by the following recursions 

(a) nt = H[i:tHt + Rt, 

(b) Kt = Ft^tHt, 

(c) yt = H[xt, (2) 

(d) xt+i = {Ft - KtH[) xt + KtVt, 

(e) Si+i = FtlltFl^ - Kt^t^K + GtQtG't, 



with starting values 

(/) xi = i?(xi) = 0, 

(g) Si = E(xix;) = l^i, 
where e^ = y^ — y^ is the y^-residuals with covariance matrix Vlt, T^t = E [(x^ — X() (x^ — X() ] is 
interpreted as the covariance matrix of the one-step state prediction errors, and Kt = E (xt+i^'t) is 
known as the Kalman gain. The notation A> means that the matrix A is nonnegative definite. 

Recursion (2e) based on the starting equation (2g) will be called periodic Riccati difference 
equation (PRDE) because in the limit, i.e. when T.t+sk converges as /c ^ oo for alH G {1, ..., S*}, 
the S'-periodic limiting solution Pt = lim T.t+sk will satisfy the following discrete-time matrix 

A;— »oo 

periodic Riccati equation (DPRE) 

Pt+^ = FtPtF[ - FtPtHt (H^PtHt + Rt)-' H[PtF[ + GtQtG[, t G {1, ..., 5}, 

which has been extensively studied (see for example Bittanti et al, 1988 for some theoretical 
aspects and Hench and Laub, 1994 for a numerical resolution). As is well known, the resolution 
of (2e) requires O(r^) operations per iteration which is computationally expensive. Furthermore, 
the solution S^ must be nonnegative definite, a property that is not easy to preserve in a numerical 
resolution of (2e). The following section proposes some recursions that avoid these drawbacks and 
may have further advantages over the Kalman filter (2). 

IL Periodic Chandrasekhar-type algorithms 

The recursions proposed in this section and which are aimed to generalize Morf et a/'s (1974) 
algorithms to the periodic case will be called analogously periodic Chandrasekhar-type equations. 
This, of course, will not mean that there is an analog of our recursions in the periodic continuous- 
time case. The derivation of our recursions is similar to its classical counterpart and is based on 
the factorization result given below (see Theorem 3.1). 



Let AgSf = S(+5 — S( denote the S'-lagged increment of the Riccati variable, for given 
Si, S2, ..., S5 > 0. Then, one can proves the following result. 
Theorem 3.1 The S-lagged increment AsT^t satisfies the following dijference equations 

As^t+i = [Ft - Kt+s^;_lsH't) [As^t + AgEtHtfi;' H^As^t] {Ft - Kt^s^-l^E'^ , (3) 

A5St+i = [Ft - KtQ;'H',) [AsJ^t - As^tHtQ^^sH'^As^t] {Ft - KtQ^'H',)' . (4) 

Proof 

i) Proof of (3) 

From (2a) we have 

^t+s = H^AsT^tHt + H^EtHt + Rt 
= H',As^tHt + nt. 

Hence 

nt+s = ^t + As^t, (5) 

where As^t = H'tAsT^tHt. Moreover, from (2e) it follows that 

Si+1+5 = Ft^t+sFt — Kt+s^t+sK^_^_g + GtQtG^, 
S,+i = FtT.tFl-KtntK[ + GtQtG[, 

where Kt = KtVl^^. Therefore, 

AsT.t+1 = FtAs^tF; - Kt+s^t+sK[+s + Kt^tK',. (6) 



On the other hand, the Kalman gain Kt may be written in a backward recursive form as follows 



Kt = {Ft^t+sHt-Ft^s^tHt)n- 



-1 



Kt+s^t+s-Ft^s^tHAn 



1-1 



Kt+s H.Ht+sHt + Rt)- FtAs^tH, 



Q 



-1 



Kt+s H.AsJ^tHt + H^EtHt + Rt - FtAs^tHt 



Q 



-1 



Kt+s^t + Kt+sH,As^tHt - FtAs^tHt 
Kt+s - (Ft - Kt+sH',) AsJ^tHtQ;'. 



n;' 



Whence 



Kt = Kt+s - AsKt, 



(7) 



with AsKt =^ (^Ft - Kt+sH[) AsT^tHt^iK 

Now replacing the latter expression of Kt in the last term of the right hand side of (6) while 
using (5), we obtain 



A.S 



s^t+i 



= FtAsT.tF[ - Kt+s^t+sK't^s + Kt^tK'^ 

= FtAsJ^tF; - Kt+s (fit + H[As^tHt) K[^s + [Kt+s - AsKt) Qt [Kt+s - AsKt 

= FtAs^tF; - Kt+s^tK[^s - Kt+sH[AsT.tHtK[^s 

+Kt+s^tK[+s - Kt+s^tAsK[ - AsKtQtK'.^s + AsKtQtAsK', 
= FtAsJ^tF; - Kt+sH[AsT.tHtK[^s " Kt+sH[As^t [Ft - Kt+sH\ 

[Ft - Kt+sH[) Asi:tHtn~^H[AsT.t [Ft - Kt+sH\ 



+ [Ft 



- [Ft - Kt+sH[) As^tHtK[^s 
= [Ft - Kt+sH'^ AsT^tFl - [Ft - Kt+sH'^ AsT.tHtK[^s 

+ [Ft - Kt+sH[) AsT.tHtn-^H[AsT.t [Ft - Kt+sH[^ 
= [Ft - Kt+sH[) [AsT.t + AsT.tHtn~^H[As^t] (Ft - Kt+sH[ 



proving (3). 



ii) Proof of (4) 

A similar argument may be used to prove (4). It suffice to express Kt-is witli respect of Kt in 
a forward recursive form as follows 



K, 



t+s 






that is 



FtJ^tHt + Ft^s^tHt) n-ls 
Kt^t + FtAs^tHt^ n;l, 
Kt {H'jltHt + R^ + FAs^tHt 
Kt [H'jlt^sHt - H[As^tHt + Rt) + FAs^tHt 
KtQt+s - KtHlAs^tHt + F^As^^Ht 
Kt + [Ft - KtH'^ AsT^tHt^jl,, 

Kt+s = Kt + AsKt, 






"t+5 



(8) 



where AsKt = [Ft - KtH[) As^tHtn;^s- 

Then, replacing the expression of Kt+s given by (8) in the second term of the right hand side 
of (6), it follows that 



A.S 



s^t+i 



FtAsT.tF[ - Kt+s^t+sK't^s + Kt^tK 

FtAsT^tFl - (Kt + AsK^ ^t+s (Kt + AsKt 

+Kt (ilt+s - H[AsT.tH^ K't 

FtAsT^tFl - Kt^t+sK't - Kt^t+sAsK't - AgKt^t+sK't - AsKt^t+sAsK't 

^Kt^t+sK't - KtH'tAsT.tHtK't 

FtAsT.tF[ - Kt^t+sAsK't - AsKtQt+sK't - AsKtQt+sAsK't - KtH'tAs^tHtK'f 



Finally, using again (8) we can write A^S^+i as follows 

As^t+i = FtAsJ^tP; - Kt ((f* - kX^ As^tHt) - (f* - kX'^ As^tHtK', 

-{{Ft- KtH'i^ AsT^tH^ n;l, {{Pt - KtH'i^ AsT^tH^ - KtH[AsT.tHtK[ 
= FtAsT^tFl - KtH[AsT.t {Ft - KtH'^ - {Pt - KtH'^ As^tHtK', 

- {Ft - KtH'i) AsT.tHtn;lsH'tAsT.t {Ft - KtH['^ - KtH[AsT.tHtK[ 
= FtAslltFl - KtH[AsT.tF[ + KtH^AsT^t {KtH[^ - FtAsT.tHtK[ 

- {Ft - KtH'i^ As^tHtn;^sHt^s^t {Ft - KtH[^ 
= {Ft - KtH'i) As^tFt - {Ft - KtH't) As^t {KtH'^ 

- {Ft - KtH'i^ As^tHtn;^sH'tAs^t {Ft - KtH'^ 

= {Ft - KtH'^ As^t {Ft - KtH'^ - {Ft - KtH',^ As^tHt^^^^sH'tAs^t {Ft - KtH[^ , 

showing (4). ■ 

Theorem 3.1 shows that AgSj may be factorized as follows 

As^t = YtMtYl, (9) 

where Mt is a square symmetric matrix, non necessarily nonnegative definite, of dimension 
rank (A5S1), which is at least equal to rank (A^S^). Indeed, from (3) we have 

rank (A^S^+i) < rank (A^S^) < ... < rank (A5S1) < r. 

This can be exploited to derive some recursions with a best computational complexity than the 
filter (2). 

Let us remark that Theorem 3.1 is not surprising since one can always write a periodically time 
varying state-space model (1) as a time-invariant state space model (see Meyer and Burrus, 1975) 
to which it may be possible to apply the standard Chandrasekhar type factorization due to Morf 
et al (1974). Nevertheless, because of the requiring increasing bookkeeping (the obtained time 



invariant system is of dimension multiplied by S) the development of a proper theory for periodic 
state-space models would be fruitful. 

Thanks to the factorization result given by Theorem 3.1, the matrices Yt and Mt given by (9) 
can be obtained recursively. The following algorithm shows that the periodic Riccati difference 
equation (2e) may be replaced by a set of recursions on Vlt, Kt, Yt and Mt with a reduction in 
computational efforts, especially when the state dimension r is much larger than m, the dimension 

of yt. 

Algorithm 3.1 The Kalman filter (2) can be replaced by a set of recursive equations containing 

(2c) and {2d) and the following recursions 

(a) Qt+s = a + mYtMtY.'Ht, 

(b) Kt+s = {Kt + FtYtMtYlHt) , 

(c) Yt+, = [Ft - Kt+s^~lsH't) Yt, 



(10) 



{d) Mt+i = Mt + MtY^Htnt'H^YtM, 



with starting values 



(e) ns = K^sHs, s=l,...,S, 

if) Ks = FsEsHs, s = l,...,S, 
where T.^,! < s < S is determined from (2e) and {2g), while Yi and Mi are obtained by factorizing 
nonuniquely 



AsSi = Fs^sF's - KsQs K's + GsQsG's - Si, 



(10^) 



as 



YiMiy;. 



Derivation (10a) is just (6) when using (9), while (106) follows from (9) and the relation 



Kt+s = {FtEtHt + FtAs^tHt 



10 



On the other hand, from (3) which we rewrite while using (9) we obtain 

As^t+i = [Ft - Kt+sH[) [^s^t + ^s^tHtn-^H[^s^i\ (f* - Kt+sH[) 

= [Ft - Kt+sH[) \YtMtYl + YtMtYlHtn-^H[YtMtYi\ [Pt - Kt+sH[ 

= {Ft - Kt+s^;lsH't) Yt [Mt + MtY;Htn-^H[YtMt) F/ [Ft - Kt+s^;l,H[)' 

= ^t+i^t+i^f+i- 

By simple identification we get (10c) and (lOd). ■ 

Note that the PRDE (2e) must be executed for 1 < s < S* to start recursions (10). However, for 
t > 5* the recursive calculation of S^ is not dealt with by the above algorithm but can be deduced 
from it through the following equation 

T^ks+s = Ss + ^ Yjs+sMjs+sY'js+s^ l< s< S. 

j=0 

Similarly to the time-invariant case (Morf et al, 1974), other forms of Algorithm 3.1 can be derived 
from Theorem 3.1. The following variant is particularly well adapted when Mi < 0, in which case 
we have Mt < for any t. This case is encountered whenever the periodic state-space model (1) 
is periodically stationary (causal) as we can see below. 

Algorithm 3.2 The following set of recursions in which (10a), (106) and (10e)-(105') (3.8a) are 
unchanged while (10c) and (lOrf) are replaced by 

(a) y,+i = {Ft - Ktn^'H^) Yt, 

(11) 

(b) Mt+i = Mt- MtY^Htn;l,HiYtMt, 

provides the same results as Algorithm 3.1. 

Derivation The derivation is similar to that of Algorithm 3.1, but is based on the factorization (4) 

rather than (3). ■ 

It is still possible to derive other forms similarly to the standard time-invariant case. The ho- 
mogenous periodic Riccati difference equation (lOrf) can be linearized using the matrix inversion 



11 

lemma (Morf et al, 1974) through which, we obtain a recursion on Mf^ rather than on Mt as 
follows 

M-\ = Mi' - YlH,n-lsH[Y,. 

It is worth noting that the periodic Chandrasekhar recursions given by Algorithm 3.1 and Algorithm 
3.2 will be preferred to the Kalman filter (2) whenever the dimension of Yt and/or Mt are signifi- 
cantly less than that of S^. These dimensions are conditioned on the good choice of the factorization 
A5S1 = YiMiY( in the initialization step which will be studied in the following section. 

III. The initialization problem 

As is well known, the most important step in the development of a Chandrasekhar algorithm 
is the initialization step because it modulates the computational complexity and hence the lack of 
numerical advantage over the Kalman filter. In our periodic case, this step depends on the relation 
between the period S, the output dimension m, and the state dimension r. First of all, suppose the 
process {xi} given by (1) is periodically stationary, that is, all the eingenvalues of the monodromy 
matrix ni=o ^s-i ^e less than unity in modulus. Let us consider two cases. 
i) Case where Sm < r: 

As pointed out in (lOg) the start up values li and Mi are determined by factorizing A5S1 as 
YiMiY(. Iterating (lO^f) 5" times as follows 

A5S1 = FsJ^sF's-KsnsK's + GsQsG's-^i 

= Fs iFs-iT^s-iF's-i — Ks-i^s-i^s-i + Gs-iQs-iG'g_i] F'g 

-KsQsK's + GsQsG's - Si 
= ^5^5-1^5-1 (FsFs-i) — FsKs-i^s-iKg^iFg + FsGs-iQs-iG's_iFg 

-KsQsK's + GsQsG's - Si 



12 



FsFs-1 iFs-2^S-2Fs^2 ~ F^S-2^S~2Kg_2 + Gs-2Qs-2Gs_2) (FsFs-l) 

-FsKs-iQs-iK's^iF's + FsGs^iQs-iG's^iF's - KsQsK's + GsQsG's - ^i 

FsFs-iFs-2^s-2 {FsFs-iFs^2) — FsFs-iKs-2^s-2 iFsFs-iKs-2 

+ {FsFs-i) Gs^2Qs-2Gg_2 (FsFs-i) 

-FsKs^iQs-i [Ks-iFs) + FsGs-iQs^iG's^iF's - KsQsK's + GsQsG's - Si 



5-1 /k-1 \ /fc-1 

k=0 \j=0 J \j=0 



^S-1 \ /5-1 

d=o / \j=o 



5-1 /fc-1 \ /fc-1 \ ' 

+ E n ^^-^ Gs^kQs-kG's_, n ^^-^ - ^1' (12) 

fc=0 \j=o / \i=o / 

and invoking the fact that under the periodic stationarity assumption, Si satisfies the following 
discrete-time periodic Lyapunov equation (DPLE) (e.g. Bittanti et al, 1988; Varga, 1997) 

(5-1 \ /5-1 \ ' 5-1 /fc-1 \ /fc-1 \ ' 

n ^s-, j^Au^^-A + E n ^^-^ Gs^kQs-kG's^k n ^^-^ > 
i=o / Vi=o / fc=o \i=o / \j=o J 

we conclude that the sum of the last three terms of the right hand-side of (12) is zero. 



Whence 



A5S1 



5-1 /fc-1 \ /fc-1 

E n ^s-, Ks-u^s-kK's^u n ^s-, 

fc=o \i=o / \i=o 

5-1 /fc-1 \ /fc-1 

fc=0 \i=0 / \i=0 



( nt 



-L 







\ 



■■■ ^rv 



L' = y,m,y;, 



(13) 



where L is given by 



L 



5-1 



Ks, FsKs-u FsFs.,Ks-2, H Fs^jK, 



13 



Clearly, with such a factorization the dimension of Mi (and hence of Mt for every t) is equal to 
mS which is less than r, the dimension of the Riccati matrix associated with the Kalman filter (2). 
Indeed, when Sm is fairly less than r, the nonhomogeneous PRDE (2e) may be replaced by the 
homogenous PRDE (life) which is of lower dimension. For instance, for m = 1, the complexity 
of solving (lOd) or (116) when using (13) as an initialization step is of order 0{Sr'^) which is 
computationally simple to solve compared to the PRDE (2e). It is still possible to improve the 
computation of (13) by alleviating the formation of the sums of products in L by using the periodic 
Schur decomposition (Bojanczyk et al, 1992; Hench and Laub, 1994). 
a) Case where Sm > r: 

In this case the latter factorization given by (13) would be inefficient since the dimension of Mj 
is greater than that of S^. Thus we have to search for another factorization. We have 

El = E (xi — Xi) (xi — Xi) = E (xix-^^j 

= E {Fsxo + Gswo) {FsXq + GsWo) 

= FsE (^xqx'qJ F's + GsE (wqIv'^j C'g 

= FsWoF's + GsQsG's- 



Therefore, 



A5S1 = FiMiFi 

= FsJ^sF;, - KsQsK's + GsQsG's - FsWoF's - GsQsG's 

= Fs{T.s~Wo)F's-KsnsK's 

= Fs (E5 - Wo) F's - {FsEsHsns') Qs (F^Esi/s^^^)' 



F. 



J:s-Wo- i^sHs) ^s i^sH, 



F's- 



14 

This allows to identify Yi and Mi as follows 

(14) 
Ml = lis -Wo- i^sHs) ^s' i^sHs)' . 

With such an initialization, the PRDE (11&) has the same dimension as that of the PRDE (2e), 
and it seems that there is no reduction in the computational cost compared to the Kalman filter. 
However, the difference from (2e) is that, unlike the S^, the Mt is not required to be nonnegative- 
definite. This helps alleviate the computational complexity of (lOd) and then (lib). 

In the matter of illustration we propose the following example which shows the impact of a good 
choice of a starting factorization on the Chandrasekhar algorithm complexity. 
Example 4.1 Consider a periodic autoregression of order 5 and period S {PARs{5)), which is 
given by the following stochastic difference equation 

Vt - (p?yt-i - 4>fyt~2 - 4*Vt-3 - (pfyt-i - 4>fyt-b = £*, (15) 

where {et} is a periodic white noise with S'-periodic variance and where the parameters (()■', 
j = 1, ..., 5 are periodic with respect of t with S. 

Setting Xi = {yt,yt-i.:,yt-4y, £t = (e*, 0,0, 0,0)' and H' = (1,0,0,0,0)', model (15) may be 

written in the state-space form 

Xi = FiXt_i + St 

(16) 

yt = H% 

so that identifying it with model (1), the dimensions r, m and d are respectively equal to 5, 1 and 
5. 

When applying the Kalman filter to model (16), the corresponding periodic Riccati equation (2e) 
is of dimension 5 (the dimension of S^) for any value of S. However, the dimension of the Riccati 
equation corresponding to the periodic Chandrasekhar filter (dimension of Mt) depends upon S. 

Let us consider two cases for S. 
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i) Case where S = 2. 

We are in the case where Sm < r. According to formula (13), we have 



A5E1 = -L 



( 1 \ 

^Z^ 



^\^ 



L' = YiM.Yl, 



J 



with L = [K2,F2Ki], Ki = FJ^iHi and ^t = H[T.tHt, t = 1,2. So, we can take Mi = 

^ ft ^ 

, from which the corresponding Riccati equation is of dimension 2, clearly lower 

than the dimension of the Riccati equation of the Kalman filter. Whence in this case the periodic 
Chandrasekhar filter is highly superior to its homologue, the Kalman one. 
ii) Case where S = 12. 

In this case, the previous factorization is inefficient since the dimension of the Chandrasekhar 
Riccati would be equal to 12, much larger than 5, the dimension of the Kalman Riccati. Nevertheless, 
we are in the case Sm > r, and according to (14), A5S1 may be factorized as YiMiY(, where 

Fi = F12 and Ml = S12 - Wq - ^i2Hi2n:,^ H[^T.[^, 

so that the Riccati equation associated with the Chandrasekhar filter has the same dimension as that 
of the Riccati equation of the Kalman filter. Moreover, the matrix Mi is not necessarily nonnegative 
definite in contrast with Si, and from this viewpoint the Chandrasekhar filter is still more suitable. 

IV. Conclusion 

In this paper the discrete-time Chandrasekhar recursions have been generalized to the periodic 
time-varying state-space case through several forms. These recursions allow in a large range of cases 
to solve the periodic Riccati difference equation with a considerable reduction in the computational 
complexity. Along similar lines to the standard time-invariant case (Morf and Kailath, 1975), a 
square root version of these recursions may be easily derived in order to improve the numerical 
stability of the proposed algorithms. Useful applications for time series analysis as well as for the 
periodic system theory can be given, in particular, we mention the likelihood evaluation of periodic 
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VARMA (Aknouche and Hamdi, 2007), the calculation of exact Fisher information matrix for 
PARMA models and the development of fast RLS algorithms for periodic systems (Bentarzi and 
Aknouche, 2006). 
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